Disordered Hubbard Model with Attraction: 
Coupling Energy of Cooper Pairs in Small Clusters 
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We generalize the Cooper problem to the case of many interacting particles in the vicinity of the 
Fermi level in the presence of disorder. On the basis of this approach we study numerically the 
variation of the pair coupling energy in small clusters as a function of disorder. We show that the 
Cooper pair energy is strongly enhanced by disorder, which at the same time leads to the localization 
of pairs. 

PACS numbers: 72.15.Rn, 71.30+h, 74.20.-z 



I. INTRODUCTION 

The superconductor-insulator transition (SIT) in dis- 
ordered films of metals has attracted widespread interest 
in recent yearsEl. The SIT, driven by adjusting some tun- 
ing parameter such as film thickness or magnetic field 
strength, is particularly interesting in two dimensions 
(2D), where both superconductivity and metallic behav- 
ior are marginal. 

In the composite bosons model, Cooper pairs are 
treated as point-like charge 2e bosons. In this picture, the 
superconducting state displays a quantum phase transi- 
tion to an insulating state characterized by a quenching 
of the condensate of composite bosonscl. In this scenario, 
the SIT is caused by the loss of phase coherence between 
the pairs in different parts of the sample, while the magni- 
tude of the pairing gap remains finite. Numerical studies 
support this scenario^. 

However, a few relevant experimental aspects of the 
SITd seem to be beyond the scope of the composite 
fermions theory. It is therefore highly desirable to study 
models where the fermionic nature of charge carriers 
is not rubbed out from the beginning. The Quantum 
Monte Carlo studies of the disordered attractive Hubbard 
model in two dimensions have supported the possibility 
of a disorder driven superconductor to insulator quantum 
phase transition.]. At the same time the mean field ap- 
proach within the Bogoliubov-de Gennes framework has 
shown that also space fluctuations of the pairing ampli- 
tude should be taken into account in order to give a full 
picture of the SITH 

In parallel a growing interest has been devoted to the 
question of what is the coupling energy of pairs placed 
in small superconducting grains, with the average ISMel 
spacing of the same order as the superconducting gapQTl 
Also the pair properties in small size samples may be re- 
lated to their properties in the localized phase, where the 
pair motion is bounded inside the localization domain. 

In this paper we study numerically the properties of 



Cooper pairs in small two-dimensional clusters with dis- 
order. We take the Hubbard attractive interaction be- 
tween fermionic particles with spin 1/2 which move in 
a two-dimensional Anderson lattice. Following the ap- 
proach introduced by Coopero, we consider some part 
of the particles below the Fermi sea as frozen, while the 
remaining particles, in the direct vicinity of the Fermi 
level, can move and interact in the presence of disor- 
der. Recently, such generalized Cooper problem has been 
consideredtlJ for the case of two particles in a disordered 
potential. Here, we further develop this approach for the 
case of many interacting Cooper pairs, that allows us to 
study the case of finite particle density. 

Our numerical studies allow us to determine the depen- 
dence of the Cooper pair coupling energy on the strength 
of disorder. They show that this coupling can be strongly 
increased by disorder, which however leads to localization 
of pairs. In the regime of weak disorder the pairs are delo- 
calized but their coupling energy is significantly reduced 
compared to the localized phase. 

The paper is organized as follows: In Section |l| we 
introduce the attractive Hubbard model and discuss the 
numerical method used to study the case of a finite parti- 
cle density. In Section III we determine the Cooper pair 



coupling energy (pairing gap) and investigate its depen- 
dence on the strength of the disorder. In Section IV 



we 



study the disorder-induced pair localization and compare 
the results with the case of noninteracting particles. In 
Section ^ we study the behavior of the superconducting 
order parameter, obtained from the pairing correlation 
function. The conclusions are presented in Section VI. 



II. MODEL AND NUMERICAL METHOD 

We study a disordered square lattice with N fermions 
on L? sites. The Hamiltonian is defined by 

H=-V ^2 c l* c i° + U ^2 n it n ih (!) 

<ij><7 itr i 
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where c\ a (cj CT ) creates (destroys) an electron at site i with 
spin a, Ti\fj — c\ 'c\ a is the corresponding occupation num- 
ber, the hopping term V between nearest neighbors lat- 
tice sites characterizes the kinetic energy, the site energies 
£i are taken from a box distribution over [-W/2, W/2], U 
measures the strength of the Hubbard attraction (U < 0) , 
and periodic boundary conditions are taken in both di- 
rections. We restrict our numerical investigations to the 
subspace with S z = for even TV (N/2 spins up and N/2 
spins down) and S z = 1/2 for odd N. 

The model (Q) at U = reduces to the one body An- 
derson model, giving localized states in two dimensions 
at the thermodynamic limitiia. At W = one gets the 
clean attractive Hubbard model, which in 2D shows a 
Kosterlitz-Thouless transition to a superconducting state 
with power-law decay of the pairing correlationald. 

We study numerically the model ([!]) for a finite density 
of interacting quasiparticles above the frozen Fermi sea: 

(i) Single particle eigenvalues E a and eigenstates (or- 
bitals) 4> a (i) (a = 1,...,L 2 ) at U = are obtained via 
numerical diagonalization of the Anderson Hamiltonian. 

(ii) The Hamiltonian ([!]) is written in the orbital basis: 

with cft aa . = <p a (i)c icr , and transition matrix elements 

Qlp = Y,4><*®fo®W)Mi)- (3) 



with E g (U) many body ground state for an attractive 
interaction U . 




FIG. 1. Dependence of the TV-body pairing energy E v on 
the number N of electrons, at different disorder strengths 
W, from W = (bottom) to W/V = 15 (top) in steps of 
AW/V — 1. The linear system size is L — 10 and the cut-off 
orbital M = 12. Here and in the following figures, U = —41/ 
and data are averaged over Nr — 100 disorder realizations. 



(iii) The Fermi sea is introduced by restricting the sums 
in (Q) to orbitals with energies above the Fermi energy 
E mp : 01,0,^,5 > mp. We consider a filling factor 
v = nip/L 2 = 1/4 (corresponding to 2mp frozen elec- 
trons due to spin degeneracy) and a finite density of N 
interacting quasiparticles above the Fermi level. 

(iv) The Slater determinants basis, built from the single 
particle orbitals cj> a , is energetically cut-off by means of 
the condition X^i( TO » — m F) _ M, with m, orbital in- 
dex for the i-th quasiparticle (m, > mp). Such a rule 
gives an effective phonon frequency u/p, oc M / L 2 . 

(v) The ground state of this tmiicated Hamiltonian is 
found via the Lanczos algorithrrU. 

In our numerical simulations we considered up to N = 
8 interacting quasiparticles in up to M = 20 orbitals. We 
checked that results are qualitatively similar under vari- 
ation of the cut-off orbital M. In the following Sections 
we present data for U = —AV, averaged over Nr = 100 
disorder realizations. 
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III. PAIRING GAP 

In order to compute the pairing energy, we first com- 
pute the total iV-body pairing energy as 

E p (N) = E g (U = 0)-E g (U), (4) 



FIG. 2. Pairing gap A{N) = E P (N) - EJN - 1) vs. dis- 
order strength W, with E p taken from Fig. ID, = 2 (stars), 
N = 4 (circles), TV = 6 (squares), and TV = 8 (diamonds). 

In Fig. |l| we show the TV-body pairing E P (N) as a 
function of the number N of interacting quasiparticles 
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above the Fermi sea, at different disorder strengths. This 
figure shows a clear even-odd effect, with a much larger 
increase of the pairing energy when TV is pair. This fact 
has a clear meaning: for TV even, it is possible to build 
a new pair, reducing the ground state energy due to the 
negative coupling U < 0. For TV odd, in an ideal BCS 
superconductor the additional particle cannot be paired 
and remains as a quasiparticle excitation. However, a 
small ground state energy reduction is still present in our 
numerical simulations, since the unpaired particle weakly 
interacts with the superconducting pairs. 

The jump in the pairing energy from odd to even num- 
ber of particles, 

A(N)=E P (N)-E P (N-1), (5) 

with TV even, can therefore be interpreted as the energy 
necessary to break a superconducting pair; in the ideal 
BCS case, this would give the superconducting energy 
gap. We note that the superconducting gap is extracted 
in a similar way in experiments with single Cooper pair 
tunneling inside superconducting islandst£IE3. 

In Fig. H we show the pairing gap A(TV) as a function 
of the disorder strength W, for TV = 2,4,6,8. We see 
that, with the exception of the first jump (TV = 2), the 
other jumps are rather similar. It is clear that A grows 
significantly with the disorder strength W. We attribute 
this effect to the fact that at strong disorder particles are 
trapped in the deepest minima of the random potential. 
Therefore the pair size becomes smaller that enhances 
the interaction between coupled particles, hence A. 

A/V 2 [ ' ' \ 




FIG. 3. Dependence of pairing gap A (average of A (TV = 2) 
and A (TV = 4)) on the disorder strength W, at L = 6 (circles), 
L = 8 (squares), and L = 10 (diamonds), keeping the ratio 
M/L 2 » 0.2. 

In Fig. H we show the dependence of A on the sys- 
tem size 6 < L < 10. Since the Debye cut-off frequency 



should be independent of the system size, we keep con- 
stant the ratio M/L 2 m 0.2. The coupling energy A be- 
comes independent of the system size at large W, while 
it is not yet saturated at small W. As the pair size is 
determined by 1/A, this means that at small W the size 
of the pair becomes comparable with the system size. 

IV. PAIR LOCALIZATION 

In order to study the localization properties of the sys- 
tem, we consider the fraction £ of the sample occupied 
by the TV-body wavefunction \^ g ): 

TV 2 

where 

Pi* = (* s |n iCT |* s ) (7) 

is the charge density of the ground state at the site i. 
With the definition (|), TV/2L 2 < £ < 1, the lower limit 
corresponding to pairs localized in a single site, the upper 
limit to complete charge derealization. 




FIG. 4. Fraction of occupied sites £ as a function of the dis- 
order strength W, for M/L 2 « 0.2, U = -4V (full symbols) 
and U = (empty symbols). Main figure: TV = 6 particles, at 
L — 6 (circles), L = 8 (squares), and L = 10 (diamonds). In- 
set: constant particle density of mobile fermions N/L 2 « 0.06; 
TV = 4, L = 8 (circles) and TV = 6, L = 10 (squares). 

In Fig.[| we show the fraction of occupied sites £ as a 
function of disorder, at different system sizes 6 < L < 1 0, 
for TV = 6 particles and a fixed Debye frequency (M / L 2 ks 
0.2). This figure gives a clear indication of the presence 
of two regimes: at small disorder the wavefunction fills 
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a large fraction of the sample (superconducting regime), 
while at large disorder £ decreases with the system size 
(£ oc l/L 2 , localized regime). In the inset of Fig. [| we 
show the parameter £ at different system sizes L = 8, 10, 
and for a constant electronic density of mobile fermions 
N/L 2 « 0.06. The two curves puts on top of each other, 
suggesting the existence of a size-independent function 
£(W) in the thermodynamic limit. The drop of £(W) 
with W demonstrates that disorder gives localization of 
Cooper pairs. 

Finally we examine the question of to what extent the 
wavefunction localization is a many body effect instead 
of a single particle Anderson localization phenomenon. 
Therefore in Fig. [| we also show the parameter £ 
at U = 0. The comparison between the interacting 
(U = —4) and the noninteracting (U — 0) case sug- 
gests that the interaction makes localization stronger, in 
agreement with results for taio particles in a three di- 
mensional random potential!!!!. This effect can he ex- 
plained qualitatively with the following argumentllll: at- 
tractive interaction creates pairs of total mass m p twice 
the electronic mass. This halves the effective hopping 
term V c g oc l/m p , thus doubling the ratio W/V e ff. We 
remark that this rough argument fails in the delocalized 
regime at small disorder, where the tendency seems to 
be reversed. ActU|al|Ly, due to localization of single par- 
ticles states in 2DLL3, for a given number N of particles 
one should get £({/ = 0) — > when L — > oo. 
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FIG. 5. Dependence of the order parameter A op on the 
number N of particles, at different disorder strengths W, from 
W = (top) to W/V = 15 (bottom) in steps of AW/V = 1. 
The linear system size is L = 10 and the cut-off orbital 
M = 12. 



V. PAIR CORRELATION 

The superconducting state can be characterized by the 
s-wave pair correlation function, 

P s (r) = (* 9 |-l^A i+r At|*A, (8) 

i 

where 

A| = (9) 

creates a pair at site i. For an s-wave superconducting 
state, 

A op = y/P s (r=(L/2,L/2)) (10) 

is the order parameter of the superconductor-insulator 
transitionB. 
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FIG. 6. Dependence of the slope a of the order parameter 
A op on the disorder strength IF; a is extracted from a linear 
fit of A op (iV) vs N, with data taken from Fig. |[ 

In Fig. H we show the dependence of the order pa- 
rameter A op on the number of interacting quasiparticles 
above the Fermi sea. The order parameter is strongly 
suppressed by disorder (see also Ref.a), an effect which 
becomes more evident with the addition of particles. We 
remark that A op shows an approximate linear increase 
with the number of pairs. This is quite natural if the 
many-.body ground state wavefunction is in the BCS 
formO-3, built from single particle eigenfunctions including 
disorderEa: 

|* s > ocTJ(l+ ST 6t)|0), (11) 

7 
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with 



dLdl 



- .* ji* 1 5 and g 1 variational parameters. Using 

the relation c ia = J2 a ^aW^ao-j after lattice and disor- 
der averaging the dominant contributions in P s (r) is pro- 
portional to the number of pairs, A op (7V) oc a% (here 
we have taken into account that g 7 in the BCS .theory 
changes smoothly with a around the Fermi levelc-3 and 
here we considered N <C 2np). Here a is a parameter 
which determines the slope of A op variation with AT. The 
order parameter decreases when an unpaired particle is 
added. In our opinion, this is due to the fact that this 
extra electron weakly interacts with the paired particles, 
reducing the pair correlation function. 

In FigJ^ we show the slope a of the linear fit of the 
order parameter A op as a function of the number N of 
quasiparticles. The suppression of this quantity with dis- 
order is evident, indicating a rather sharp crossover from 
a superconducting to an insulating behavior in our finite 
size lattice. 



VI. CONCLUSIONS 

In this paper we have investigated the localization of 
Cooper pairs for small clusters in a two-dimensional dis- 
ordered substrate. We have shown that the Cooper pair 
coupling energy displays an even-odd asymmetry: this 
parity effect survives also in the presence of disorder. 
The pairing gap is strongly enhanced by disorder, which 
at the same time leads to localization of Cooper pairs 
(gapped insulator). Therefore, in the insulating regime, 
the breaking of Cooper pairs should enhance transport. 
This is consistent with the resistivity drop .observed in 
experiments with an applied magnetic fielcOtj, which 
might signal the crossover from a Cooper pair insulator 
to an electronic insulator. 
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